Exploratory analysis

Load data

vcs <- rio::import("data_complete_all_anonym_allZ_unfiltered.xlsx")
var_label(vcs$dominance) <- "Dominance"
var_label(vcs$neuro) <- "Neuroticism"
var_label(vcs$agree) <- "Agreeableness"
var_label(vcs$extra) <- "Extraversion"
var_label(vcs$openn) <- "Openness"
var_label(vcs$consc) <- "Conscientiousness"
var_label(vcs$soi_full) <- "Unrestricted sociosexuality"
var_label(vcs$f0) <- "Voice pitch"
var_label(vcs$pf) <- "Formants"
vcs$sex_c <- vcs$sex
vcs$sex <- factor(if_else(vcs$sex == 1, "male", "female"))
contrasts(vcs$sex) <- contr.helmert(2)
var_label(vcs$sex) <- "Sex"
set.seed(1)
var_label(vcs$age) <- "Age"
vcs <- vcs %>% 
  mutate(
    age = if_else(dataset == 9, 20, age)/10,
    age_se = if_else(dataset == 9, 3, 0.5)/10
  )
warmup <- 2000
iter <- warmup + 2000
chains <- 4
control <- list(adapt_delta = 0.99)
priors <- c(
  prior(normal(0, 3), class = b)
)

variable_labels <- c("Intercept"= "Intercept",  "sex1" = "Sex [male]", "bsp_meageage_se" = "Age±SE", "f0" = "Voice pitch (f0)", "pf" = "Formant position (Pf)", "age" = "Age")
effect_labels <- c("b_Intercept"= "Intercept",  "b_sex1" = "Sex [male]", "bsp_meageage_se" = "Age±SE", "b_f0" = "Voice pitch (f0)", "b_pf" = "Formant position (Pf)", "b_age" = "Age")

Missingness patterns

codebook::md_pattern(vcs %>% select(dominance, extra, behavior, sex, age, f0, pf, dataset))

## # A tibble: 6 x 9
##   description                       age    pf    f0 behavior extra dominance var_miss n_miss
##   <chr>                           <dbl> <dbl> <dbl>    <dbl> <dbl>     <dbl>    <dbl>  <dbl>
## 1 Missing values per variable         4     4     7      242   803      1256     2316   2316
## 2 Missing values in 0 variables       1     1     1        1     1         1        0    969
## 3 Missing values in 2 variables       1     1     1        1     0         0        2    760
## 4 Missing values in 1 variables       1     1     1        1     1         0        1    267
## 5 Missing values in 2 variables       1     1     1        0     1         0        2    196
## 6 9 other, less frequent patterns     7     7     5        2     3         2       28     49

vcs <- vcs %>% filter(!is.na(f0), !is.na(pf), !is.na(age))

Neuroticism

Visually diagnose non-linearity

ggplot(vcs, aes(f0, neuro)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

ggplot(vcs, aes(pf, neuro)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

ggplot(vcs, aes(age*10, neuro)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  scale_x_continuous("Age") +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

Decisions upon visual diagnosis:

Models

h0 <- brm(neuro ~ sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/neuro/h0") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

h1_simple <- brm(neuro ~ f0 + pf + sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/neuro/h1_simple") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

compare_models <- loo_compare(h0, h1_simple)
compare_models

##           elpd_diff se_diff
## h1_simple  0.0       0.0   
## h0        -0.7       2.4

sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)

	Neuroticism
Predictors	Estimates	CI (95%)
Intercept	0.07	-0.25;0.36
Voice pitch (f0)	0.17	0.03;0.32
Formant position (Pf)	-0.05	-0.17;0.07
Sex [male]	-0.12	-0.29;0.05
Age±SE	-0.02	-0.14;0.09
Random Effects
σ²	0.01
τ₀₀	0.99
ICC	0.01
N _dataset	7
Observations	1434
Marginal R² / Conditional R²	0.088 / 0.088

Model details

summary(h1_simple)

## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be careful when analysing the results!
## We recommend running more iterations and/or setting stronger priors.

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: neuro ~ f0 + pf + sex + me(age, age_se) + (1 | dataset) 
##    Data: vcs (Number of observations: 1434) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.14      0.08     0.03     0.33 1.00     2758     3241
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       0.07      0.15    -0.25     0.36 1.00     4674     5172
## f0              0.17      0.08     0.03     0.32 1.00     6520     5752
## pf             -0.05      0.06    -0.17     0.07 1.00    10242     5972
## sex1           -0.12      0.09    -0.29     0.05 1.00     5816     5934
## meageage_se    -0.02      0.06    -0.14     0.09 1.00     4782     5651
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.96      0.02     0.92     0.99 1.00    15341     5954
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

conditional_effects(h1_simple)

# conditional_smooths(h1_after_diagnosis)
equiv_test <- equivalence_test(h1_simple, range = c(-0.1, 0.1))
equiv_test_table <- equiv_test %>% 
  mutate(Parameter = recode(Parameter, !!!effect_labels)) %>% 
  mutate(
    ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
    HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>% 
  select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")

Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).
Parameter	H0	inside ROPE	89% HDI
Intercept	Undecided	49%	[-0.19; 0.30]
Voice pitch (f0)	Undecided	15%	[ 0.05; 0.28]
Formant position (Pf)	Undecided	81%	[-0.15; 0.05]
Sex [male]	Undecided	39%	[-0.26; 0.01]
Age±SE	Undecided	95%	[-0.12; 0.06]

plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)

## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).

Openness

Visually diagnose non-linearity

ggplot(vcs, aes(f0, openn)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

ggplot(vcs, aes(pf, openn)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

ggplot(vcs, aes(age*10, openn)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  scale_x_continuous("Age") +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

Decisions upon visual diagnosis:

Models

h0 <- brm(openn ~ sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/openn/h0") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

h1_simple <- brm(openn ~ f0 + pf + sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/openn/h1_simple") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

compare_models <- loo_compare(h0, h1_simple)
compare_models

##           elpd_diff se_diff
## h1_simple  0.0       0.0   
## h0        -0.9       2.4

sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)

	Openness
Predictors	Estimates	CI (95%)
Intercept	-0.15	-0.47;0.18
Voice pitch (f0)	-0.19	-0.34;-0.03
Formant position (Pf)	-0.01	-0.15;0.12
Sex [male]	-0.26	-0.45;-0.07
Age	0.04	-0.07;0.15
Random Effects
σ²	0.01
τ₀₀	0.99
ICC	0.01
N _dataset	7
Observations	1434
Marginal R² / Conditional R²	0.024 / 0.024

Model details

summary(h1_simple)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: openn ~ f0 + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 1434) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.20      0.11     0.05     0.48 1.00     1458     2308
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.15      0.16    -0.47     0.18 1.00     4555     4251
## f0           -0.19      0.08    -0.34    -0.03 1.00     4960     5149
## pf           -0.01      0.07    -0.15     0.12 1.00     3767     4639
## sex1         -0.26      0.09    -0.45    -0.07 1.00     3502     4478
## age           0.04      0.06    -0.07     0.15 1.00     6261     6161
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.99      0.02     0.95     1.02 1.00     7539     5573
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

conditional_effects(h1_simple)

# conditional_smooths(h1_after_diagnosis)
equiv_test <- equivalence_test(h1_simple, range = c(-0.1, 0.1))
equiv_test_table <- equiv_test %>% 
  mutate(Parameter = recode(Parameter, !!!effect_labels)) %>% 
  mutate(
    ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
    HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>% 
  select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")

Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).
Parameter	H0	inside ROPE	89% HDI
Intercept	Undecided	33%	[-0.42; 0.09]
Voice pitch (f0)	Undecided	9%	[-0.31;-0.06]
Formant position (Pf)	Undecided	94%	[-0.13; 0.10]
Sex [male]	Rejected	0%	[-0.41;-0.11]
Age	Undecided	89%	[-0.04; 0.13]

plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)

## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).

Conscientiousness

Visually diagnose non-linearity

ggplot(vcs, aes(f0, consc)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

ggplot(vcs, aes(pf, consc)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

ggplot(vcs, aes(age*10, consc)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  scale_x_continuous("Age") +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

Decisions upon visual diagnosis:

Models

h0 <- brm(consc ~ sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/consc/h0") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

h1_simple <- brm(consc ~ f0 + pf + sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/consc/h1_simple") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

compare_models <- loo_compare(h0, h1_simple)
compare_models

##           elpd_diff se_diff
## h0         0.0       0.0   
## h1_simple -1.9       0.5

sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)

	Conscientiousness
Predictors	Estimates	CI (95%)
Intercept	-0.14	-0.57;0.32
Voice pitch (f0)	-0.03	-0.17;0.12
Formant position (Pf)	0.02	-0.11;0.15
Sex [male]	-0.12	-0.29;0.06
Age	0.00	-0.11;0.11
Random Effects
σ²	0.11
τ₀₀	0.90
ICC	0.10
N _dataset	7
Observations	1434
Marginal R² / Conditional R²	0.119 / 0.119

Model details

summary(h1_simple)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: consc ~ f0 + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 1434) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 7) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.44      0.18     0.22     0.89 1.00     1881     3332
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.14      0.22    -0.57     0.32 1.00     2714     3640
## f0           -0.03      0.07    -0.17     0.12 1.00     4862     5011
## pf            0.02      0.07    -0.11     0.15 1.00     5811     5587
## sex1         -0.12      0.09    -0.29     0.06 1.00     4439     4693
## age           0.00      0.06    -0.11     0.11 1.00     7718     5262
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.94      0.02     0.91     0.98 1.00     7783     5907
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

conditional_effects(h1_simple)

# conditional_smooths(h1_after_diagnosis)
equiv_test <- equivalence_test(h1_simple, range = c(-0.1, 0.1))
equiv_test_table <- equiv_test %>% 
  mutate(Parameter = recode(Parameter, !!!effect_labels)) %>% 
  mutate(
    ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
    HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>% 
  select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")

Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).
Parameter	H0	inside ROPE	89% HDI
Intercept	Undecided	31%	[-0.47; 0.24]
Voice pitch (f0)	Undecided	87%	[-0.15; 0.09]
Formant position (Pf)	Undecided	94%	[-0.09; 0.12]
Sex [male]	Undecided	40%	[-0.27; 0.02]
Age	Accepted	100%	[-0.09; 0.09]

plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)

## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).

SOI Desire

Visually diagnose non-linearity

ggplot(vcs, aes(f0, desire)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

ggplot(vcs, aes(pf, desire)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

ggplot(vcs, aes(age*10, desire)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  scale_x_continuous("Age") +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

Decisions upon visual diagnosis: Visual diagnostic show no clear sign of nonlinearities. We will fit linear effects for f0, pf, and age. It is not necessary to model measurement error in age, as dataset 9 (where age was not noted individually) is not included in these analyses. Visually, there could be an interaction between sex and f0.

Models

h0 <- brm(desire ~ sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/desire/h0") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

h1_simple <- brm(desire ~ f0 + pf + sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/desire/h1_simple") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

compare_models <- loo_compare(h0, h1_simple)
compare_models

##           elpd_diff se_diff
## h1_simple  0.0       0.0   
## h0        -5.4       3.6

sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)

	desire
Predictors	Estimates	CI (95%)
Intercept	0.35	0.03;0.69
Voice pitch (f0)	-0.22	-0.34;-0.11
Formant position (Pf)	0.07	-0.05;0.20
Sex [male]	0.26	0.12;0.40
Age	-0.13	-0.22;-0.04
Random Effects
σ²	-0.01
τ₀₀	1.01
ICC	-0.01
N _dataset	9
Observations	2000
Marginal R² / Conditional R²	0.164 / 0.164

Model details

summary(h1_simple)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: desire ~ f0 + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 2000) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 9) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.38      0.13     0.21     0.69 1.00     2120     3184
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.35      0.17     0.03     0.69 1.00     2670     3457
## f0           -0.22      0.06    -0.34    -0.11 1.00     5673     5127
## pf            0.07      0.06    -0.05     0.20 1.00     5522     5205
## sex1          0.26      0.07     0.12     0.40 1.00     4950     5207
## age          -0.13      0.05    -0.22    -0.04 1.00     7094     5433
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.92      0.01     0.89     0.95 1.00     7229     4298
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

conditional_effects(h1_simple)

# conditional_smooths(h1_after_diagnosis)
equiv_test <- equivalence_test(h1_simple, range = c(-0.1, 0.1))
equiv_test_table <- equiv_test %>% 
  mutate(Parameter = recode(Parameter, !!!effect_labels)) %>% 
  mutate(
    ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
    HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>% 
  select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")

Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).
Parameter	H0	inside ROPE	89% HDI
Intercept	Undecided	0%	[ 0.09; 0.62]
Voice pitch (f0)	Rejected	0%	[-0.32;-0.13]
Formant position (Pf)	Undecided	67%	[-0.03; 0.18]
Sex [male]	Rejected	0%	[ 0.14; 0.37]
Age	Undecided	25%	[-0.20;-0.05]

plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)

## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).

SOI Attitude

Visually diagnose non-linearity

ggplot(vcs, aes(f0, attitude)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

ggplot(vcs, aes(pf, attitude)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

ggplot(vcs, aes(age*10, attitude)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  scale_x_continuous("Age") +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

Models

h0 <- brm(attitude ~ sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/attitude/h0") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

h1_simple <- brm(attitude ~ f0 + pf + sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/attitude/h1_simple") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

compare_models <- loo_compare(h0, h1_simple)
compare_models

##           elpd_diff se_diff
## h1_simple  0.0       0.0   
## h0        -7.1       4.4

sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)

	attitude
Predictors	Estimates	CI (95%)
Intercept	0.02	-0.34;0.39
Voice pitch (f0)	-0.24	-0.36;-0.12
Formant position (Pf)	-0.07	-0.20;0.06
Sex [male]	-0.07	-0.21;0.07
Age	-0.00	-0.09;0.09
Random Effects
σ²	0.08
τ₀₀	0.92
ICC	0.08
N _dataset	9
Observations	1999
Marginal R² / Conditional R²	0.132 / 0.132

Model details

summary(h1_simple)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: attitude ~ f0 + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 1999) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 9) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.43      0.15     0.25     0.80 1.00     2186     3484
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.02      0.19    -0.34     0.39 1.00     2477     3149
## f0           -0.24      0.06    -0.36    -0.12 1.00     5375     5107
## pf           -0.07      0.07    -0.20     0.06 1.00     6077     6108
## sex1         -0.07      0.07    -0.21     0.07 1.00     4715     4801
## age          -0.00      0.05    -0.09     0.09 1.00     7425     5033
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.93      0.01     0.91     0.96 1.00     7581     5266
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

conditional_effects(h1_simple)

# conditional_smooths(h1_after_diagnosis)
equiv_test <- equivalence_test(h1_simple, range = c(-0.1, 0.1))
equiv_test_table <- equiv_test %>% 
  mutate(Parameter = recode(Parameter, !!!effect_labels)) %>% 
  mutate(
    ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
    HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>% 
  select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")

Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).
Parameter	H0	inside ROPE	89% HDI
Intercept	Undecided	48%	[-0.27; 0.31]
Voice pitch (f0)	Rejected	0%	[-0.34;-0.15]
Formant position (Pf)	Undecided	71%	[-0.17; 0.04]
Sex [male]	Undecided	68%	[-0.19; 0.05]
Age	Accepted	100%	[-0.08; 0.07]

plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)

## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).

Sociosexuality

Visually diagnose non-linearity

ggplot(vcs, aes(f0, soi_full)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

ggplot(vcs, aes(pf, soi_full)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

ggplot(vcs, aes(age*10, soi_full)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  scale_x_continuous("Age") +
  geom_smooth(method = "gam", colour = "black", formula = y ~ s(x)) +
  geom_smooth(aes(colour = sex, fill = sex), method = "gam",
              formula = y ~ s(x)) +
  scale_color_viridis_d("Sex") +
  scale_fill_viridis_d("Sex")

Models

h0 <- brm(soi_full ~ sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/soi_full/h0") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

h1_simple <- brm(soi_full ~ f0 + pf + sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/soi_full/h1_simple") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

compare_models <- loo_compare(h0, h1_simple)
compare_models

##           elpd_diff se_diff
## h1_simple   0.0       0.0  
## h0        -10.0       5.0

sjPlot::tab_model(h1_simple, bpe = "mean",show.ngroups = T,ci.hyphen = ";",
                  pred.labels = variable_labels)

	Unrestricted sociosexuality
Predictors	Estimates	CI (95%)
Intercept	-0.20	-0.60;0.19
Voice pitch (f0)	-0.29	-0.41;-0.17
Formant position (Pf)	-0.01	-0.14;0.12
Sex [male]	-0.03	-0.17;0.12
Age	0.09	-0.01;0.18
Random Effects
σ²	0.07
τ₀₀	0.93
ICC	0.07
N _dataset	9
Observations	1995
Marginal R² / Conditional R²	0.160 / 0.160

Model details

summary(h1_simple)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: soi_full ~ f0 + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 1995) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 9) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.46      0.15     0.27     0.85 1.00     1910     3239
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.20      0.20    -0.60     0.19 1.00     2156     3365
## f0           -0.29      0.06    -0.41    -0.17 1.00     4808     4766
## pf           -0.01      0.06    -0.14     0.12 1.00     6259     5195
## sex1         -0.03      0.07    -0.17     0.12 1.00     4258     4948
## age           0.09      0.05    -0.01     0.18 1.00     7240     5169
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.92      0.01     0.89     0.95 1.00     7124     4467
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

conditional_effects(h1_simple)

# conditional_smooths(h1_after_diagnosis)
equiv_test <- equivalence_test(h1_simple, range = c(-0.1, 0.1))
equiv_test_table <- equiv_test %>% 
  mutate(Parameter = recode(Parameter, !!!effect_labels)) %>% 
  mutate(
    ROPE_Percentage = sprintf("%2.f%%", 100*ROPE_Percentage),
    HDI = sprintf("[% .2f;% .2f]", HDI_low, HDI_high)) %>% 
  select(Parameter, H0 = ROPE_Equivalence, `inside ROPE` = ROPE_Percentage, `89% HDI` = HDI)
kable(equiv_test_table, caption = "Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).")

Test whether the 89% highest-density interval (HDI) is fully within the region of practical equivalence (ROPE) of [-0.1;0.1] that we set (where all continuous variables, except age are standardised).
Parameter	H0	inside ROPE	89% HDI
Intercept	Undecided	26%	[-0.51; 0.13]
Voice pitch (f0)	Rejected	0%	[-0.39;-0.19]
Formant position (Pf)	Undecided	96%	[-0.12; 0.09]
Sex [male]	Undecided	89%	[-0.14; 0.09]
Age	Undecided	60%	[ 0.02; 0.16]

plot(equiv_test) + scale_y_discrete("Effect", breaks = names(effect_labels), labels = effect_labels)

## Warning: Removed 3516 rows containing non-finite values (stat_density_ridges).

Extraversion, adjusted for dominance

extra_adj_for_dominance <- brm(extra ~ dominance + f0 + pf + sex + me(age, age_se) + (1 | dataset), data = vcs,
          iter = iter + 1000, warmup = warmup + 1000, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/extra/extra_adj_for_dominance2") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2")

## Warning: Rows containing NAs were excluded from the model.

## Warning in system2(CXX, args = ARGS): error in running command

## Warning in file.remove(c(unprocessed, processed)): cannot remove file '/var/folders/tr/
## lgwjmdgd72bckr1hj5gd_7640000gq/T//RtmpP8sN6o/filed8fc6940eed1.stan', reason 'No such file or directory'

## Warning in system2(CXX, args = ARGS): error in running command

## Warning in file.remove(c(unprocessed, processed)): cannot remove file '/var/folders/tr/
## lgwjmdgd72bckr1hj5gd_7640000gq/T//RtmpP8sN6o/filed8fc33e44653.stan', reason 'No such file or directory'

## Warning: The largest R-hat is 1.1, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#r-hat

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess

## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#tail-ess

extra_adj_for_dominance

## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be careful when analysing the results!
## We recommend running more iterations and/or setting stronger priors.

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: extra ~ dominance + f0 + pf + sex + me(age, age_se) + (1 | dataset) 
##    Data: vcs (Number of observations: 970) 
## Samples: 4 chains, each with iter = 5000; warmup = 3000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.25      0.24     0.03     0.91 1.00     1401     1844
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       0.27      0.22    -0.12     0.70 1.00     2210     2299
## dominance       0.49      0.03     0.44     0.54 1.00     8896     6242
## f0             -0.08      0.08    -0.24     0.08 1.00     4535     4679
## pf              0.17      0.08     0.02     0.32 1.00     5156     5544
## sex1           -0.08      0.09    -0.27     0.10 1.00     3758     4835
## meageage_se    -0.09      0.05    -0.20     0.01 1.00     4368     4818
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.79      0.02     0.76     0.82 1.00     9652     5775
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

dominance_adj_for_extra <- brm(dominance ~ extra  + f0 + pf + sex + age + (1 | dataset), data = vcs,
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/dominance/dominance_adj_for_extra") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 
dominance_adj_for_extra

## Warning: There were 3 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: dominance ~ extra + f0 + pf + sex + age + (1 | dataset) 
##    Data: vcs (Number of observations: 970) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.51      0.41     0.16     1.66 1.00     2008     1880
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.22      0.35    -0.97     0.44 1.00     2004     1655
## extra         0.54      0.03     0.49     0.60 1.00     9622     5817
## f0           -0.14      0.08    -0.30     0.02 1.00     5687     5227
## pf           -0.11      0.08    -0.27     0.04 1.00     7063     5231
## sex1         -0.17      0.10    -0.37     0.02 1.00     5034     4957
## age           0.08      0.05    -0.02     0.19 1.00     8144     5164
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.83      0.02     0.79     0.87 1.00    10065     4955
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Extraversion, separate by sex

h1_men <- brm(extra ~ f0 + pf + me(age, age_se) + (1 | dataset), data = vcs %>% filter(sex == "male"),
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/extra/h1_men") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

h1_women <- brm(extra ~ f0 + pf + me(age, age_se) + (1 | dataset), data = vcs %>% filter(sex != "male"),
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/extra/h1_women") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

h1_men

## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be careful when analysing the results!
## We recommend running more iterations and/or setting stronger priors.

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: extra ~ f0 + pf + me(age, age_se) + (1 | dataset) 
##    Data: vcs %>% filter(sex == "male") (Number of observations: 617) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 5) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.34      0.26     0.04     0.99 1.00     1598     2017
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept      -0.12      0.32    -0.72     0.52 1.00     4357     4488
## f0             -0.21      0.13    -0.47     0.06 1.00    10846     5489
## pf              0.04      0.10    -0.15     0.24 1.00    10328     5930
## meageage_se    -0.09      0.08    -0.25     0.07 1.00     5646     5778
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.00      0.03     0.94     1.06 1.00    13493     6418
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

h1_women

## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be careful when analysing the results!
## We recommend running more iterations and/or setting stronger priors.

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: extra ~ f0 + pf + me(age, age_se) + (1 | dataset) 
##    Data: vcs %>% filter(sex != "male") (Number of observations: 816) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 5) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.50      0.30     0.19     1.27 1.00     2079     2745
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       0.39      0.34    -0.29     1.07 1.00     2544     3368
## f0             -0.24      0.10    -0.43    -0.05 1.00     9072     5621
## pf              0.14      0.09    -0.05     0.32 1.00     8805     5928
## meageage_se    -0.11      0.08    -0.26     0.05 1.00     6496     5416
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.96      0.02     0.91     1.00 1.00     9974     5586
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Agreeableness, separate by sex

h1_men <- brm(agree ~ f0 + pf + me(age, age_se) + (1 | dataset), data = vcs %>% filter(sex == "male"),
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/agree/h1_men") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

h1_women <- brm(agree ~ f0 + pf + me(age, age_se) + (1 | dataset), data = vcs %>% filter(sex != "male"),
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/agree/h1_women") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

h1_men

## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be careful when analysing the results!
## We recommend running more iterations and/or setting stronger priors.

## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: agree ~ f0 + pf + me(age, age_se) + (1 | dataset) 
##    Data: vcs %>% filter(sex == "male") (Number of observations: 617) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 5) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.47      0.30     0.16     1.24 1.00     2048     2830
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       0.08      0.34    -0.57     0.77 1.00     2760     3319
## f0              0.22      0.12    -0.02     0.46 1.00    11523     5334
## pf             -0.04      0.10    -0.23     0.15 1.00    10635     5748
## meageage_se    -0.02      0.08    -0.17     0.14 1.00     4997     5371
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.92      0.03     0.87     0.98 1.00    11879     6533
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

h1_women

## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be careful when analysing the results!
## We recommend running more iterations and/or setting stronger priors.

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: agree ~ f0 + pf + me(age, age_se) + (1 | dataset) 
##    Data: vcs %>% filter(sex != "male") (Number of observations: 817) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 5) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.74      0.39     0.32     1.78 1.00     2403     3966
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept       0.24      0.42    -0.60     1.08 1.00     2754     3158
## f0             -0.11      0.10    -0.30     0.09 1.00     9247     5453
## pf              0.09      0.09    -0.09     0.27 1.00    10273     5582
## meageage_se    -0.11      0.08    -0.27     0.05 1.00     8799     5832
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.95      0.02     0.90     1.00 1.00    11155     5888
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Dominance, separate by sex

h1_men <- brm(dominance ~ f0 + pf + age + (1 | dataset), data = vcs %>% filter(sex == "male"),
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/dominance/h1_men") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

h1_women <- brm(dominance ~ f0 + pf + age + (1 | dataset), data = vcs %>% filter(sex != "male"),
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/dominance/h1_women") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

h1_men

## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: dominance ~ f0 + pf + age + (1 | dataset) 
##    Data: vcs %>% filter(sex == "male") (Number of observations: 494) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 3) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.80      0.75     0.18     2.92 1.00     1452     1720
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.39      0.69    -1.70     0.85 1.00     2233     1950
## f0           -0.26      0.16    -0.56     0.05 1.00     7084     5782
## pf            0.06      0.19    -0.32     0.44 1.00     7318     5055
## age           0.04      0.09    -0.13     0.21 1.00     6316     4626
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.97      0.03     0.91     1.04 1.00     7379     5184
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

h1_women

## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: dominance ~ f0 + pf + age + (1 | dataset) 
##    Data: vcs %>% filter(sex != "male") (Number of observations: 491) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 3) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.68      0.63     0.13     2.33 1.00     1887     2941
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.22      0.55    -0.93     1.27 1.00     2874     2798
## f0           -0.28      0.12    -0.52    -0.03 1.00     6677     5433
## pf           -0.05      0.11    -0.27     0.16 1.00     6240     4875
## age           0.03      0.09    -0.16     0.21 1.00     5275     5369
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.97      0.03     0.91     1.04 1.00     6843     5405
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

SOI behavior, separate by sex

h1_men <- brm(behavior ~ f0 + pf + age + (1 | dataset), data = vcs %>% filter(sex == "male"),
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/behavior/h1_men") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

h1_women <- brm(behavior ~ f0 + pf + age + (1 | dataset), data = vcs %>% filter(sex != "male"),
          iter = iter, warmup = warmup, chains = chains, cores = chains,
          prior = priors,
          control = control, save_mevars = TRUE, 
          file = "models/behavior/h1_women") %>% 
  add_criterion("loo") %>% 
  add_criterion("bayes_R2") %>% 
  add_criterion("loo_R2") 

h1_men

## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.99 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: behavior ~ f0 + pf + age + (1 | dataset) 
##    Data: vcs %>% filter(sex == "male") (Number of observations: 790) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 6) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.47      0.23     0.21     1.07 1.00     1947     2574
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.96      0.31    -1.54    -0.36 1.00     3287     3753
## f0           -0.29      0.12    -0.51    -0.06 1.00     6649     5601
## pf            0.05      0.15    -0.24     0.35 1.00     5791     5444
## age           0.25      0.07     0.12     0.37 1.00     6570     5551
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.93      0.02     0.88     0.98 1.00     7034     4764
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

h1_women

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: behavior ~ f0 + pf + age + (1 | dataset) 
##    Data: vcs %>% filter(sex != "male") (Number of observations: 1206) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~dataset (Number of levels: 6) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.43      0.23     0.20     0.98 1.00     1853     2329
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.93      0.27    -1.45    -0.39 1.00     2948     3625
## f0           -0.15      0.07    -0.29    -0.01 1.00     7158     5076
## pf           -0.05      0.07    -0.19     0.09 1.00     6429     5018
## age           0.44      0.07     0.30     0.57 1.00     6851     5510
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.92      0.02     0.88     0.96 1.00     7620     5324
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Summary of results

Models

openn <- readRDS("models/openn/h1_simple.rds")
neuro <- readRDS("models/neuro/h1_simple.rds")
consc <- readRDS("models/consc/h1_simple.rds")
desire <- readRDS("models/desire/h1_simple.rds")
attitude <- readRDS("models/attitude/h1_simple.rds")
soi_full <- readRDS("models/soi_full/h1_simple.rds")
theme_set(theme_cowplot())

modify <- function(outcome) {
  list(scale_y_discrete(outcome, breaks = names(effect_labels), labels = effect_labels),
  guides(fill = "none"),
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0, 0.5, 0, 0.5), "cm"))
)
}
remove_x_axis <- list(
  scale_x_continuous("", limits = c(-0.5, 0.4)),
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_blank(),
        axis.line.x = element_blank())
)

p_extra <- plot(equivalence_test(neuro, parameters = c("b_f0", "b_pf", "b_sex1"))) + modify("Neuroticism") + remove_x_axis
p_openn <- plot(equivalence_test(openn, parameters = c("b_f0", "b_pf", "b_sex1"))) +
  modify("Openness") + remove_x_axis
p_consc <- plot(equivalence_test(consc, parameters = c("b_f0", "b_pf", "b_sex1")))  +
  modify("Conscientiousness") + remove_x_axis
p_attitude <- plot(equivalence_test(attitude, parameters = c("b_f0", "b_pf", "b_sex1"))) +
  modify("SOI attitude") + remove_x_axis
p_desire <- plot(equivalence_test(desire, parameters = c("b_f0", "b_pf", "b_sex1"))) +
  modify("SOI desire") + remove_x_axis
p_soi_full <- plot(equivalence_test(soi_full, parameters = c("b_f0", "b_pf", "b_sex1"))) +
  modify("Sociosexuality") + xlim(c(-0.5, 0.4))
legend <- get_legend(
  p_extra + 
    guides(fill = guide_legend(nrow = 1)) +
    theme(legend.position = "bottom", legend.justification = "center")
)

## Warning: Removed 2637 rows containing non-finite values (stat_density_ridges).

neg_spacer <- -0.155
plot_grid(ncol = 1,  rel_heights = c(1, neg_spacer, 1, neg_spacer, 1,neg_spacer, 1,neg_spacer, 1, neg_spacer, 1.1, .15), align = "v", axis = "b",
          p_extra, NULL, p_openn,NULL, p_consc,NULL, p_attitude, NULL, p_desire, NULL, p_soi_full, legend)

## Warning: Removed 2637 rows containing non-finite values (stat_density_ridges).

## Warning: Removed 2637 rows containing non-finite values (stat_density_ridges).

## Warning: Removed 2637 rows containing non-finite values (stat_density_ridges).

## Warning: Removed 2637 rows containing non-finite values (stat_density_ridges).

## Warning: Removed 2637 rows containing non-finite values (stat_density_ridges).

## Warning: Removed 2637 rows containing non-finite values (stat_density_ridges).

ggsave("exploratory_rope.pdf", width = 8, height = 10)

Tabular

get_coefs <- function(model) {
model_summary <- summary(model)
nonvarying <- model_summary$fixed %>% 
  as.data.frame() %>% 
  rownames_to_column("term") %>% 
  mutate(group = "non-varying",
         outcome = as.character(model_summary$formula[[1]][[2]]))
varying <- model_summary$random %>% 
  map(~ rownames_to_column(as.data.frame(.), "term")) %>% 
  bind_rows(.id = "group") %>% 
  left_join(model_summary$ngrps %>% 
    as_tibble() %>% 
    gather(group, n), by = "group") %>% 
  mutate(#group = paste0(group, " (n=", n, ")"),
         outcome = as.character(model_summary$formula[[1]][[2]]))

two_digits <- function(x) { sprintf("%.2f", x) }
coefs <- bind_rows(nonvarying, 
                   varying)
coefs <- coefs %>% select(group, outcome, term, Estimate, `l-95% CI`, `u-95% CI`)
coefs$Estimate = two_digits(coefs$Estimate)
coefs$`l-95% CI` = two_digits(coefs$`l-95% CI`)
coefs$`u-95% CI` = two_digits(coefs$`u-95% CI`)
coefs
}

coefs <- bind_rows(get_coefs(neuro), 
                   get_coefs(openn), get_coefs(consc), get_coefs(desire),
                   get_coefs(attitude),get_coefs(soi_full))

## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be careful when analysing the results!
## We recommend running more iterations and/or setting stronger priors.

wide_coefs <- coefs %>% 
  mutate(Estimate = paste0(Estimate, " [", `l-95% CI`, ";", `u-95% CI`, "]")) %>% 
    select(group, outcome, term, Estimate) %>% 
  mutate(term = recode(term, "meageage_se" = "Age(±SE)", "age" = "Age(±SE)"),
         outcome = fct_inorder(outcome),
         term = fct_inorder(term)) %>% 
  spread(outcome, Estimate) %>% 
  arrange(desc(group), term) %>% 
  mutate(group = fct_inorder(group))

outcome_names <- c("Neuroticism" = "neuro", "Conscientiousness" = "consc",
                   "Openness" = "openn", "Sociosexuality" = "soi_full",
                   "SOI desire" = "desire", "SOI attitude" = "attitude")

library(kableExtra)
wide_coefs %>% 
  mutate(term = str_replace_all(term, variable_labels)) %>% 
  rename(!!!syms(outcome_names), Term = term) %>% 
  select(-group) %>% 
  kable(caption = "Exploratory analysis results") %>%
  pack_rows(index = table(wide_coefs$group)) %>% 
  column_spec(column = 1, width = "3.3cm") %>% 
  column_spec(column = 2:(ncol(wide_coefs)-1), width = "3cm") %>% 
  add_header_above(c(" " = 1, "Estimated effect on each outcome [95% CI]" = ncol(wide_coefs)-2)) %>% 
  footnote("Estimated associations between outcomes, for which we had made no predictions.")

Exploratory analysis results
	Estimated effect on each outcome [95% CI]
Term	Neuroticism	Openness	Conscientiousness	SOI desire	SOI attitude	Sociosexuality
non-varying
Intercept	0.07 [-0.25;0.36]	-0.15 [-0.47;0.18]	-0.14 [-0.57;0.32]	0.35 [0.03;0.69]	0.02 [-0.34;0.39]	-0.20 [-0.60;0.19]
Voice pitch (f0)	0.17 [0.03;0.32]	-0.19 [-0.34;-0.03]	-0.03 [-0.17;0.12]	-0.22 [-0.34;-0.11]	-0.24 [-0.36;-0.12]	-0.29 [-0.41;-0.17]
Formant position (Pf)	-0.05 [-0.17;0.07]	-0.01 [-0.15;0.12]	0.02 [-0.11;0.15]	0.07 [-0.05;0.20]	-0.07 [-0.20;0.06]	-0.01 [-0.14;0.12]
Sex [male]	-0.12 [-0.29;0.05]	-0.26 [-0.45;-0.07]	-0.12 [-0.29;0.06]	0.26 [0.12;0.40]	-0.07 [-0.21;0.07]	-0.03 [-0.17;0.12]
Age(±SE)	-0.02 [-0.14;0.09]	0.04 [-0.07;0.15]	0.00 [-0.11;0.11]	-0.13 [-0.22;-0.04]	-0.00 [-0.09;0.09]	0.09 [-0.01;0.18]
dataset
sd(Intercept)	0.14 [0.03;0.33]	0.20 [0.05;0.48]	0.44 [0.22;0.89]	0.38 [0.21;0.69]	0.43 [0.25;0.80]	0.46 [0.27;0.85]
Note:
Estimated associations between outcomes, for which we had made no predictions.

Scatterplots

outcome_names <- c("Neuroticism" = "neuro", "Conscientiousness" = "consc",
                   "Openness" = "openn", "Sociosexuality" = "soi_full",
                   "SOI desire" = "desire", "SOI attitude" = "attitude")

neuro <- ggplot(vcs, aes(f0, neuro)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(aes(colour = sex, fill = sex), method = "lm") +
  scale_color_manual("Sex", values = c("#40A5BF", "black")) +
  scale_fill_manual("Sex", values = c("#40A5BF", "black"))  +
  guides(color = F, fill = F) +
  xlab("") +
  ylab("Neuroticism") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))

consc <- ggplot(vcs, aes(f0, consc)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(aes(colour = sex, fill = sex), method = "lm") +
  scale_color_manual("Sex", values = c("#40A5BF", "black")) +
  scale_fill_manual("Sex", values = c("#40A5BF", "black"))  +
  guides(color = F, fill = F) +
  xlab("") +
  ylab("Conscientiousness") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))

openn <- ggplot(vcs, aes(f0, openn)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(aes(colour = sex, fill = sex), method = "lm") +
  scale_color_manual("Sex", values = c("#40A5BF", "black")) +
  scale_fill_manual("Sex", values = c("#40A5BF", "black"))  +
  guides(color = F, fill = F) +
  xlab("") +
  ylab("Openness") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))


soi_full <- ggplot(vcs, aes(f0, soi_full)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(aes(colour = sex, fill = sex), method = "lm") +
  scale_color_manual("Sex", values = c("#40A5BF", "black")) +
  scale_fill_manual("Sex", values = c("#40A5BF", "black"))  +
  guides(color = F, fill = F) +
  xlab("") +
  ylab("Sociosexuality") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))


desire <- ggplot(vcs, aes(f0, desire)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(aes(colour = sex, fill = sex), method = "lm") +
  scale_color_manual("Sex", values = c("#40A5BF", "black")) +
  scale_fill_manual("Sex", values = c("#40A5BF", "black"))  +
  guides(color = F, fill = F) +
  xlab("Voice pitch (f0)") +
  ylab("SOI Desire") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))


attitude <- ggplot(vcs, aes(f0, attitude)) + 
  geom_jitter(aes(colour = sex), alpha = 0.3) +
  geom_smooth(aes(colour = sex, fill = sex), method = "lm") +
  scale_color_manual("Sex", values = c("#40A5BF", "black")) +
  scale_fill_manual("Sex", values = c("#40A5BF", "black"))  +
  guides(color = F, fill = F) +
  xlab("Voice pitch (f0)") +
  ylab("SOI Attitude") +
  theme(panel.spacing = unit(c(0), "cm"),
        plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))

plot_grid(ncol = 2,  openn, neuro, consc, soi_full, attitude, desire,
          align = "hv")

## Warning: Removed 796 rows containing non-finite values (stat_smooth).

## Warning: Removed 796 rows containing missing values (geom_point).

## Warning: Removed 796 rows containing non-finite values (stat_smooth).

## Warning: Removed 796 rows containing missing values (geom_point).

## Warning: Removed 796 rows containing non-finite values (stat_smooth).

## Warning: Removed 796 rows containing missing values (geom_point).

## Warning: Removed 235 rows containing non-finite values (stat_smooth).

## Warning: Removed 235 rows containing missing values (geom_point).

## Warning: Removed 231 rows containing non-finite values (stat_smooth).

## Warning: Removed 231 rows containing missing values (geom_point).

## Warning: Removed 230 rows containing non-finite values (stat_smooth).

## Warning: Removed 230 rows containing missing values (geom_point).

ggsave("scatterplots_exploratory.pdf", width = 8, height = 8)

Exploratory analysis

Ruben Arslan

3/5/2019

Load data

Missingness patterns

Neuroticism

Visually diagnose non-linearity

Models

Openness

Visually diagnose non-linearity

Models

Conscientiousness

Visually diagnose non-linearity

Models

SOI Desire

Visually diagnose non-linearity

Models

SOI Attitude

Visually diagnose non-linearity

Models

Sociosexuality

Visually diagnose non-linearity

Models

Extraversion, adjusted for dominance

Extraversion, separate by sex

Agreeableness, separate by sex

Dominance, separate by sex

SOI behavior, separate by sex

Summary of results

Models

Tabular

Scatterplots